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Abstract Neurons are characterised by a morphological structure unique amongst 
biological cells, the core of which is the dendritic tree. The vast number of dendritic 
geometries, combined with heterogeneous properties of the cell membrane, continue 
to challenge scientists in predicting neuronal input-output relationships, even in the 
case of sub-threshold dendritic currents. The Green's function obtained for a given 
dendritic geometry provides this functional relationship for passive or quasi-active 
dendrites and can be constructed by a sum-over-trips approach based on a path in- 
tegral formalism. In this paper, we introduce a number of efficient algorithms for 
realisation of the sum-over-trips framework and investigate the convergence of these 
algorithms on different dendritic geometries. We demonstrate that the convergence 
of the trip sampling methods strongly depends on dendritic morphology as well as 
the biophysical properties of the cell membrane. For real morphologies, the number 
of trips to guarantee a small convergence error might become very large and strongly 
affect computational efficiency. As an alternative, we introduce a highly-efficient ma- 
trix method which can be applied to arbitrary branching structures. 
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1 Introduction 

Discovered more than a century ago by Santiago Ramon y Cajal [1], dendrites form 
the vast majority of the surface area of a neuron, with the dendritic trees of some mo- 
toneurons representing up to 97% of total neuronal surface area and 75% of the total 
neuronal volume [2]. These complex branching structures are responsible for trans- 
ferring electrical activity between synapses and the soma. As technology evolved, 
interest in dendrites began to gather momentum, with the invention of sharp mi- 
cropipette electrodes in the early 1950s allowing intracellular recordings to be made. 
It was the breakthrough work of Wilfrid Rail [3] on the application of cable theory 
to dendritic modelling that provided significant insight into the role of dendrites in 
processing synaptic inputs, the historical perspective of which is summarised in a 
book by Segev, Rinzel and Shepherd [4] . Recent experimental and theoretical studies 
reinforce the fact that dendritic morphology and membrane properties play an impor- 
tant role in dendritic integration [5,6]. We refer the reader to the book Dendrites [7], 
devoted exclusively to these formations and revealing their biological complexity at 
different scales. 

It has also been known for some time that nonlinear voltage-gated ion channels 
are present in the dendrites of various types of neurons [8], and many recent dendritic 
models are constructed by combining the linear (passive) properties of dendrites to- 
gether with nonlinear (active) dynamics of membrane channels. Although the nonlin- 
ear properties of ion channels contribute considerably to neuronal input-output rela- 
tions, it is important to recognise that the passive properties of dendritic membranes 
provide the fundamental core for signal filtration and integration, and thus remain an 
essential component in understanding electrical signalling in dendrites [9] . 

When branched dendritic fibres are modelled by passive cable equations, the volt- 
age response across the branching structure for any form of applied current can be 
calculated via a convolution operation, as long as the Green's function for the given 
dendritic tree is found. This approach provides an alternative to the compartmental 
method, based on the discrete spatial approximation of the potential [10, 1 1]. It is not 
always trivial to construct such a Green's function for realistic dendritic geometries. 
Arbitrarily-branching systems are inherently difficult to solve, a fact recognised early 
by Rail, who proposed a method of mapping the branching structure onto an equiv- 
alent cylinder provided that certain geometrical restrictions were satisfied [12]. The 
work of Koch and Poggio [13], based on the graphical calculus of Butz and Cowan 
[14], focused on the calculation of the response function for complete dendritic trees 
in the Laplace (frequency) domain. Later, Rail's method of equivalent cylinders was 
extended by releasing the constraints on diameters of individual branches and by 
constructing the Green's function, again in the Laplace domain [15, 16]. An alter- 
native method for constructing the Green's function for a branching structure with 
a shunted soma was proposed by Evans and coauthors [17-19]. In this series of pa- 
pers, the response function was found in the form of an eigenfunction expansion, 
which converges particularly rapidly for large times. For smaller times, a Laplace- 
domain series solution provides better accuracy, agreeing well with an earlier "sum- 
over-trips" method for constructing the Green's function directly in the time domain, 
proposed by Abbott et al. [20]. This sum-over-trips framework is built on a path in- 
tegral formulation and enables the calculation of the Green's function on an arbitrary 
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dendritic geometry as a convergent infinite series solution. Cao and Abbott [21] pre- 
sented an algorithm for a computational realisation of the sum-over-trips approach, 
based on the division of trips into four classes. They applied this algorithm to a num- 
ber of sample dendritic trees, the largest of which had 22 branches, in contrast to real 
dendritic geometries, which might have more than 400 terminals alone [22], with a 
large variation in branch length. This complexity in neuronal morphologies across 
different types of neurons is expected to affect the convergence of computational im- 
plementations of the sum-over-trips framework. 

In this paper, we introduce and investigate a number of efficient algorithms for 
calculating the Green's function on dendritic trees using the sum-over-trips formal- 
ism. In Sect. 2, we review the theoretical framework and the four-classes algorithm 
of Cao and Abbott [21], and introduce alternative algorithms for the sum-over- trips 
method in Sect. 3. We begin with a modification on the four-classes algorithm aimed 
at improving its time complexity by developing a formal grammar to derive the trips. 
Then, a length-priority ordering of the trips using Eppstein's algorithm [23] for find- 
ing the k shortest trips on a graph is proposed. We also derive a stochastic approach 
for sampling trips on the tree based on a Monte-Carlo approach. Finally, a highly- 
efficient deterministic method for discretised tree structures is described. We assess 
the convergence of the introduced algorithms on different dendritic geometries in 
Sect. 4, where we also compare the delay and attenuation of voltage spread on four 
reconstructed dendritic morphologies. Finally, in Sect. 5, we provide a discussion of 
our results, as well as possible extensions of this work. 



2 The Sum-over-trips Framework 

We consider a dendritic branching structure with the dynamics of the membrane volt- 
age on a finite branch i described by the passive cable equation. An external current 
Ij(t) is injected at a location y on branch j. The transmembrane voltage across the 
dendritic tree is then described by the following set of equations: 

^-V U 0<x<Li,i^j (1) 

K 

^-Vj+Six-yVjit), 0<x<Lj. (2) 

K 

Here, a\ is the diameter of branch i (measured in um), R a is the specific cytoplasmic 
resistivity (in £"2 cm), C is the specific membrane capacitance (in uFcm -2 ), and R is 
the resistance across one unit area of passive membrane (in Q cm 2 ). Introducing the 
electrotonic space constant A.,- = ^/aiR/(4R a ), the membrane time constant r = RC 
and the diffusion coefficient D\ = X 2 /r, Eqs. (1) and (2) can be rewritten as 

0<x<Li,i?j, (3) 

+ — ^(x - y)Ij(t), 0<x<Lj. (4) 
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In addition to these equations, the appropriate boundary conditions must be specified 
at all branching nodes and terminals: continuity of the potential across a node and 
Kirchoff 's law of conservation of current. Continuity of the potential requires that, 
for all pairs of branches m and n attached to a node, 



V m (L m ,t) = V n (0,t), 



where the distal end of branch m is connected to the proximal end of branch n. Con- 
servation of current for the same node imposes 



E 



1 3V m 



dx 



x=L n 



x=0 



where r n = AR a /(ita^) is the axial resistance on branch n (in £"2 cm -1 ), and each sum 
is over all branches connected to this node either with their distal or proximal ends. 
At individual terminals, we can either impose a closed-end boundary condition, 



dVk 
dx 



-0, 



=L k 



or an open-end boundary condition, 



V k (L k ,t) = 0, 

where x = L& is a terminal on branch k. 

When the injected current has the form of a delta pulse, that is, Ij(t) = 8(t), the 
solution to Eqs. (3) and (4) is the Green's function Gij(x, y, t) which can be found 

as 

1 v-^ 

Gij(x, y, t) = > A trip Goo (L trip , t), (5) 

TtajC 

J trips 

where the sum is over all trips (more formally, graph-theoretic walks), starting at 
x and finishing at y, and describes the time-course of the membrane voltage at the 
location x on branch i in response to the injected current at the location y on branch 
j, where i can be taken to equal j if desired. The function takes the form 

Goo (L trip , t) = -Ve-^^/^e-^, (6) 

where L tr i p = Ltrip A*, y/^j) is me length of a trip along the tree that starts at 
point x/Xi on branch i and ends at point y/kj on branch j. Note that the length of 
each branch needs to be scaled by its own electrotonic space constant before L tr i p is 
calculated for Eq. (6). A constructed trip is allowed to reflect on or pass through any 
node on the tree an arbitrary number of times. The coefficients A tr i p depend on the 
constructed trip and are determined according to the following rules [20] : 

• From any starting point, A tr i p = 1 . 
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• For every node at which the trip passes from branch m to branch k where m ^ k, 
A tr ip is multiplied by a factor 2 . 

• For every node at which the trip reflects along on a node back onto the same branch 
n, A tr ip is multiplied by a factor 2p n — 1. 

• For every terminal, A tr ip is multiplied by + 1 for the closed-end boundary condition 
or by — 1 for the open-end boundary condition. 

When the electrical properties of the cell membrane are identical for all branches, the 
factors pk are defined as 

3/2 

^ = ^372' (7) 

where the sum is over all branches m connected to the node. When the parameters R 
and R a vary from branch to branch, the expression (7) must be modified: 

fan)- 1 _ 

/ , m V A m' m) 

However, note that the sum-over- trips method for constructing the Green's function 
in the time domain only works for uniform characteristic time constant r across the 
entirety of the dendritic tree. The generalisation of this framework to support a quasi- 
active membrane, instead of a passive membrane, releases this restriction and dif- 
ferent cell membrane properties can be chosen on each branch [24]. However, this 
means that the construction of the Green's function as an infinite series solution can 
only be performed in the Laplace domain. 

Knowing the Green's function for a given dendritic structure allows one to find 
the voltage response along the entire tree. By finding G//(x, y, t) for the ordered 
pair (x, y), the Green's function Gji (y,x,t) can be found using a simple reciprocity 
identity: 

Gji(y 9 x,t) = ^p-Gij(x 9 y,t). (9) 

Dm 

The voltage response can then be found for an arbitrary number of different discrete 
inputs as a sum of convolution integrals: 



Vi ( XJ )=\ Gij(x,Xj,t-s)Ij(s)ds, (10) 



where Xj is a location of a stimulus Ij(t) on branch j. 

The Green's function calculated by Eq. (5) for any branching structure with finite 
length branches includes an infinite number of terms. It is possible to show that this 
infinite series solution converges faster than e - ^, for sufficiently-high k, the number 
of nodes visited by the trip. We demonstrate this in the Appendix for an arbitrary tree 
with nodes of degree d = 3 or less. This generalises Abbott's convergence analysis 
[25], where it was shown that, for an infinite binary tree, the sum of coefficients A tr ip 
is Oil) for trips visiting any number of nodes. 
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Fig. 1 Four classes of trips 
between two points. Class 1 is 
the most direct trip, leaving x in 
the direction of y and not going 
past y before finishing (xBy). 
Class 2 leaves x in the other 
direction, but finishes when it 
meets y (xABy). Class 3 moves 
from x towards y, but goes past 
y and changes direction 
immediately after passing it 
before finishing (xBCy). 
Class 4 trips move from x first 
away from point y and pass y 
before reflecting on the next 
node and finishing (xABCy) 



2.1 Four-Classes Algorithm 

Cao and Abbott [21] introduced an algorithm for constructing the Green's function 
using the sum-over- trips method. Their algorithm is based on finding the shortest 
trips between any two points of measurement x and current injection y on a tree. 
Starting from the most direct, shortest trip from x to y, passing through the mini- 
mum number of nodes, four classes of trips are defined by allowing a trip to leave 
the point x in either direction and approach y from either direction along their re- 
spective branches. These initial trips, therefore, form the first and shortest trips in 
their respective classes; longer trips are generated incrementally from these. New ad- 
ditional trips can pass the points x and y any number of times and are allowed to 
change direction at any node. We will refer to this method as the four-classes algo- 
rithm. 

Figure 1 shows a model branching structure with two points x and y, and the four 
shortest classes of trips between them. Trips are represented as sequences of node 
identifiers, beginning and ending with x and y respectively. For example, we denote 
a trip from x to y via nodes A, B and C by its full description, xABCy. 

From these main trips, the four-classes algorithm generates all x -> y trips by 
inserting what are described as "excursions" into the trips. If A and B are adjacent 
nodes in a tree, then an excursion could be added to the trip x By to generate the trip 
xBABy, representing a reflection on node B towards A, reflecting at the terminal 
A back towards B, passing through this node and finally onto point y. This process 
can be iterated indefinitely, generating a trip with two more nodes each time. If this 
process is applied to every node on every trip with n and n + 1 nodes, then every trip 
with n + 2 and n + 3 nodes will be generated. Thus, from the four shortest x -> y trips 
on the tree, it is possible to construct all trips up to some threshold number of nodes 
in length explicitly. The lengths and coefficients of these trips can then be calculated 
from their full trip descriptions, allowing the Green's function given by Eq. (5) to be 
approximated. 
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3 Algorithmic Realisations 

Here, we suggest possible modifications to the four-classes algorithm of Cao and 
Abbott [21] as well as introduce novel alternative algorithms for the sum-over- trips 
formalism. 

3.1 Formal Language Theory Approach 

The four-classes algorithm generates duplicate trips [21], which must then be re- 
moved by a binary search through the list of existing trips for every new trip gener- 
ated, which takes 0(k logk) time overall, where k is the number of trips constructed. 
There are two different mechanisms by which duplicate trips are generated, and both 
mechanisms can be eliminated by applying simple restrictions to the choice of excur- 
sions applicable to a trip. As an example of the first, for the tree in Fig. 1 , it is possible 
to generate the trip xB ABC By in two different ways from the shortest Class 1 trip, 
xBy: 

Excursion Excursion 

t-» B — > BAB n a n B — >BCB rt A rt ^ rt 

xBy > xBABy > xBABCBy, 

Excursion Excursion 

t-» B — >BC B r* r* B — > B A B rt a n ^ rt 

xBy > xBCBy > xBABCBy. 

Due to the fact that the excursion may be added at any step in the trip (at the first or 
second B), the same trip may be generated multiple times. If we insist that excursions 
cannot be added at any step that precedes the excursion most recently added to the 
trip, this can be prevented. In the theory of context-free grammars, this is equivalent 
to requiring a leftmost derivation. We can represent this using a symbol to separate 
the mutable and immutable parts of the trip: 

Excursion Excursion 
B—^B\AB B^BICB 

x\By — ► xB\ABy — ► xBAB\CBy, 

Excursion 

x\By B ^ BlCB y x B\CBy / > xBABCBy. 

The second mechanism by which duplicate trips are produced is the addition of 
excursions along the same branch, starting from either end. In the structure in Fig. 1, 
we have both A —> A\BA and B —> B \ AB. Hence, in spite of the leftmost derivation 
rule, we can generate x ABA By in two different ways (brackets added for clarity): 

x\ABy -> x(A\BA)By, 
x\ABy -> xA(B\AB)y. 

This problem can be avoided by assigning each branch a direction. If the branch AB 
is given the direction B A, then the excursion A -> A \ B A is disallowed. The choice of 
direction for each branch is unambiguous on acyclic structures: apart from the branch 
on which x is found, each branch must be directed away from x . The branch upon 
which x resides is directed away from y. This ensures that each node has a sequence 
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of excursions that allow the algorithm to generate trips including it. The allocation of 
direction to each branch can be performed before the process of generating trips and 
may coincide with finding the four main classes of trips. These modifications require 
that the graph be acyclic, since "away from a point" is not generally definable on a 
graph with cycles. There do exist cyclic graphs for which an unambiguous grammar 
can generate the language of x —> y trips, but these are not relevant to the study of 
single dendritic trees. 

The two presented modifications of the four-classes algorithm are sufficient to 
prevent the generation of any duplicate trips, without any trips being missed. To- 
gether, they provide an unambiguous context-free grammar generating the language 
of x —> y trips. 

3.2 Length-Priority Method 

Since the coefficients A tr i p decay at most with e _Ltri P (although the number of trips 
increases with e Ltri ?), the dominating term in the Green's function (5) is the exponen- 

-L 2 

tial decay e tn P in Goo- The four-classes algorithm [21] does not generate trips in 
monotonic order in length, since trips are constructed by adding the same excursion 
to all four classes of trips. If, for example, a Class 2 trip is significantly longer than 
its Class 1 counterpart, due to x being along a long edge but close to a node, then a 
longer Class 2 trip will be generated before a potentially shorter Class 1 trip having an 
additional excursion on a shorter branch. In general, trips are likely to be disordered 
in length if the branches upon which x or y reside are substantially longer than at 
least one other branch on the tree, or if x or y are much closer to one of their adjacent 
nodes than to the other. 

Here, we propose to realise the sum-over-trips framework by a length-priority 
method. In this implementation, trips are generated and the corresponding terms 
^tripGoo(£trip, 0 are added to the infinite series solution (5) in monotonic order in 
length Ltrip. This is achieved by incorporating Eppstein's algorithm [23] for finding 
the k shortest trips on a graph in G(m + n log ft + k) time, with n being the number 
of nodes and m the number of edges on the branching structure. 

Both the four-classes algorithm and the improvements described in the language- 
theoretic approach rely on storing trips explicitly as sequences of nodes. This con- 
sumes 0(kn) space and time for k trips with n nodes but allows on-the-fly calculation 
of coefficients A tr i p . This is contrary to Eppstein's algorithm [23], which stores trips 
using an implicit representation and allows us to find the k shortest trips implicitly 
using only (9(1) space and time for each trip. The current implementation, based on 
Eppstein's algorithm, requires G(kn) time to calculate coefficients despite the sav- 
ings on space due to the implicit trip representation. However, Eppstein provides a 
method for computing any property that can be described by a monoid in 0(1) time 
per trip. Such a description of coefficient calculation exists, and its use would supple- 
ment the current G(kn) to G(k) decrease in space requirements with an analogous 
decrease in time complexity. The savings in space already allow the length-priority 
method to scale better than the four-classes algorithm. 
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3.3 Monte-Carlo Method 

The path integral formulation of the solution to the cable equation introduced by Ab- 
bott et al. [20] is derived via consideration of a Feynman-Kac representation of the 
solution in terms of random walkers on the dendritic geometry. Hence, it is natural to 
consider Monte-Carlo approaches to evaluating this path integral. Instead of a length- 
ordered series solution as provided by the length-priority approach, the Green's func- 
tion (5) can be constructed using a stochastic algorithm. The aim of this approach 
is to sample from trips x — >► y in such a way that the probabilistically more likely 
samples coincide with the trips that contribute most to the series solution (5). 

To motivate this Monte-Carlo approach, let us consider a linear diffusion equation 
along an infinite one-dimensional cable, 

dg d 2 g 

Yt =D- 2 , telO.Tl, (ID 

satisfying the initial condition Q(x, 0) = 8(x — y). Analogously, a diffusion process 
for the state variable X t can be defined by the stochastic equation 

dX t = V2DdW t , (12) 

with the Wiener process W t and the initial condition Xo = y. It is well known that 
Eq. (11) is the Kolmogorov equation of the diffusion process (12), that is, the time 
evolution equation of the probability density for the state of the diffusion (12). On 
the one hand, solution of (11) via classical numerical or analytical methodology in- 
forms the probability density of X t \ on the other hand, repeated sampling from (12) 
converges upon the solution g(x,t) of (\\). This method of sampling from random 
walks can also be applied for arbitrary geometries by setting the appropriate boundary 
conditions at the branching nodes and terminals. Knowing Qjj (ij,/)ona branching 
structure, we can easily find a solution of the cable equation on this geometry using 
the relation 

G i j(x,y,t) = g ij (x,y,t)e- t/T . 

Because the path integral form of the solution is equivalent to the expectation of a 
function on random walks upon the branching points of the dendritic tree, reduction 
of the random walk problem from the complete continuous space geometry of the 
neuron to the discrete topology of the branching points of the neuron gives a con- 
siderable efficiency saving to a Monte-Carlo solver. We introduce a parameter k mSiX , 
the maximum number of discrete hops on nodes for which we wish to calculate the 
expectation. The maximum number is based upon the effective maximum range of 
diffusion during the interval [0, t]. Then, we generate a realisation of a random walk 
on the nodes, 

Q) = (CD\,(L>2, •••,^ m J, (13) 

where each cok is a label identifying a particular node. For trips x — > y we select 
co\ such that it is either of the two nodes adjacent to the branch containing x, with 
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equal probability. By indexing a branch between two nodes, cok-i and co^, as the £th 
branch, subsequent steps are performed with the transition probability 

P(a)k\cok-i) = Pk, 2<&<& max , (14) 

where pk is given by (7). This connects the Monte-Carlo method to the earlier dis- 
cussed path enumeration methods. Here, we introduce two auxiliary functions, 0 and 
a, of sub walks of co. The first is a function indicating whether a sub walk of k steps 
on a realisation co is a valid trip, and is defined by 



0(y, k, t, co) = 



Goo(L(y, k,co),t), if cok-i and co\ are the nodes adjacent to y, 
0, otherwise, 



where k is the number of hops on nodes in the sub walk and L(y,k,co) is the length 
of the subwalk. The other auxiliary function a is defined as 



a(k, co) ■ 



1, if* =1,2, 

2, ifo)k-2^Q>k, 
(2pk - 1)/Pk, if cok-2 = m, 

1 , if at a closed terminal (this takes priority). 



The relevant function on paths can be defined as a composite of the auxiliary func- 
tions described above: 



A(y,co,t) = y^20(y,fc, t,co)\ 



k=l M=l 



The expectation of A with respect to the random walk (14) is equivalent to solving 
for the path integral, up to some value of & max at time t : 

E P [A(y, co, t)] = ^2 P(o>)A(y, co, t) 

CO 

kmax I k \ 

co k=l \i=l I 

^max / k \ 

= J2J2 2P(^)G 0O (L(y, k, co), t) Y\a(i, co) 



v. k=i \i=l 

at k 



^ ^ ^trip G oo (^trip > 0 > 



trips 

where P(&>) is the probability of the realisation co, and E denotes the expectation 
operator. Therefore, the Monte-Carlo strategy is to sample, sequentially or in parallel, 
the random function A in order to construct this expectation. 
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3.4 Matrix Method 

An alternative method of constructing the sum-over-trips series solution is by group- 
ing trips by their lengths: 



where the sum over / is over all possible trip lengths L tr ip. On a dendritic tree, dis- 
cretised as in compartmental models [10] or in a manner similar to the discretisation 
of the tree into segments in NEURON [26], grouping trips according to their lengths 
allows us to count the number of trips of a given length / without having to explicitly 
construct them. 

This method uses a modified directed edge adjacency matrix of the discretised 
tree in order to compute the sum of coefficients of trips of a given length. It requires 
all compartments to have the same fixed length Ax, although this restriction can be 
relaxed in a generalisation presented at the end of this section. The extremities of 
compartments define the position of nodes; there is a directed edge in both directions 
between adjacent nodes. 

We begin by defining V as the set of nodes and £ as the set of directed edges in the 
discretised tree. Edges are ordered pairs of nodes: e = (u, v) e £ is a directed edge 
from u to v, with u,v eV. For any edge e = (u,v),we denote the reverse edge by 
e' = (v, u). Trips are taken to begin from a point x along a starting edge s = (s\, S2) 
and end at a point y along a goal edge g = (g\ , #2) for s, g e £. We say that x e s 
or x e (s\, S2) if x resides along edge s = (s\,S2). Based on the locations of x e s 
and y e g, the orientations of s and g are defined such that the shortest x —> y trip 
satisfies x -> S2 —> • • • — > gi -> y. Therefore, the shortest x —> y trip always starts 
on edge s, that is, in the s\ — >► S2 direction, and approaches y along the edge g, in 
the g\ -> g2 direction. This is equivalent to a Class 1 trip; Class 2 trips leave x along 
the s' = (s2, s\) edge, arrive at y e g; Class 3 trips go from x e s to y e g'\ Class 4 
trips, finally, go from x e s f to y e g' . The locations of the points x e s and y e g 
along their respective edges are given as a fraction of the branch length such that 
x Ax denotes the distance from x to ^2 and yAx is the distance between node g\ and 
point y. We distinguish between k, the number of edges travelled in a particular trip, 
from the length of the trip Ltrip. Because x and y reside along their respective edges, 
the total length of a trip that travels along k edges is less than if the full distance along 
k edges had been travelled. That is, L tr i p < kAx for any combination of x, y and for 



The aim of the matrix method is to group all trips starting on a given edge and fin- 
ishing on a target edge by their lengths, Ltrip, and calculate the sum of the coefficients 
Atrip of those trips for each particular group instead of calculating coefficients indi- 
vidually for each trip. The sums of coefficients are computed simultaneously for trips 
ending on all edges, starting at a point x e (s\,S2), allowing us to compute Gij for 
all j, for all x on edge i and for all y on any edge, in a single run of the calculation. 




trips 



trips with 

^trip— I 



all*. 
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We define the coefficients function c s k : £ — > R, s e £, as the sum of all coefficients 
A^p which begin at point x e s and travel over k edges, finishing on a given edge g\ 



c k(s)= ^ Atri P' xes,yeg. 



trips 
in k jumps 

Because the set of edges £ is both ordered and finite, then c s k = (c s k (e\ ),..., 
c s k (e\£\)) g R'^' can be thought of as a vector, where e\ e £ for z = 1, . . . , \£ |. The 
ith element of the vector c s k corresponds to the sum of coefficients A tr i p for all trips 
originating at x on s and ending along the ith edge e\, having travelled over k edges. 
The vector c\ consists mostly of zeros, with a one only in the entry correspond- 
ing to the edge s, as the coefficient of moving in this direction remains 1, while 
all other moves are invalid by travelling over only one edge, and hence have coeffi- 
cient 0. 

We can now define a matrix Q e R^ x ^ such that 

Q h ci=c k +i. (15) 

Q is a modified form of the edge- adjacency matrix, where instead of containing ones 
to denote edge adjacency and zero otherwise, it contains the coefficient taken in mov- 
ing from one edge to another. The entries of Q can be computed based on the mor- 
phology of the graph. If the jth entry corresponds to edge (u, v) and the ith entry 
to edge (v, w), then the entry Qji is the coefficient taken when moving from branch 
(u, v) to (v, w). In the general case, these numerical values must be determined for 
each entry. However, in the simplified case where the radii on all branches are equal 
and all nodes have degree d = 1 or d = 3, the matrix Q can be constructed according 
to 



— if j = (u, v) and i = (i>, u), where v is a node of degree d = 3, 
1, if j = (u, v) and i = (v, u), where v is a closed terminal (d = 1), 
^, if j = (u, v) and i = (v, w), where u ^w, 
0, otherwise. 

(16) 



Note that the above rules apply to the transpose of Qij . 

Thus, knowing the matrix Q from the dendritic geometry and the vector c\ from 
the starting edge ^, it is possible to construct the sum of c s k (g) terms, for all k < k mSiX , 
equal to the sum of coefficients for all trips travelling up to k mSLX edges, from x e s to 
y G g. However, by considering trips moving from x in one direction only and arriv- 
ing at y from only one direction, we have calculated the coefficients of just Class 1 
trips. In order to find coefficients for the remaining three classes, we must also com- 
pute c s k (g), c s k (g') and c s k (g f ). These can be found in the same way as above. Using 
(15), the Green's function in (5) can therefore be written as 



Gij(x, y,t)= ^ A tripGoo(£trip, t) 



trips 
x to y 
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^max 

= J2[(Q k - l c\)G 00 (L 1 (k),t) 

+ (2*- 1 cf) g G 00 (L 2 (*),*) 

+ (Q k - 1 c[) gf G 00 (L 3 (k),t) 

+ (Q k - l c() gf G OQ (L 4 (k),t)l (17) 

where (Qc\) g is the gth element of the matrix- vector product of Q and c\. Lengths 
L\ , . . . , L4 are the lengths of Class 1 to Class 4 trips, respectively, and are defined as 

Li(k) = Ax(2(k- l)+x + y), 
L 2 (k) = Ax(2k -x +y), 
L 3 (k) = Ax(2k + x -y), 
U(k) = Ax(2(k + 1) - x - y). 

By selecting a small Ax , branches may be approximated by a discretisation using 
an integer number of edges of length Ax. As in compartmental models, this allows 
the full morphology of the dendritic tree to be approximated, in a trade-off between 
high speed (large Ax) and accuracy (small Ax). As Ax -> 0, however, this approach 
tends to the computational complexity of naively integrating the cable equation us- 
ing numerical methods. As in numerical simulations, where reducing Ax in order to 
increase accuracy brings about a necessary and associated change in At, the same 
is true of the matrix method: selecting a small Ax and hence increasing \£\, implies 
that £ max must be increased. 

This algorithm can be generalised to accept several discrete edge lengths 
Ax 1, . . . , Ax n , at an exponential cost in the number of different lengths n, allow- 
ing "caricature" neurons to be constructed from a small number of different edge 
lengths. Our description of this method is focused on the case where x and y are 
located on different branches. For computations where x and y are required to exist 
on the same edge, the edges can be discretised such that x and y appear on different 
segments. In all cases with bounded node degree, Q is a sparse matrix with only a 
few entries per row, and 0(|£|) entries altogether, making the complexity for the cal- 
culation of all coefficients 0(|£|&max) by using highly-efficient sparse linear algebra 
algorithms. 

3.4.1 Example calculation 

Here, we demonstrate an example realisation of the matrix method for a dendritic 
structure of three branches of equal length Ax, shown in Fig. 2. In this symmetrical 
case, the matrix Q is very small and can be constructed by hand. We place the point 
of measurement x along edge s = (A, B) and the point of current injection y along 
g = (B,D). 

We begin by ordering the edge pairs as follows: (A, B), (B, A), (B, C), (C, B), 
(B, D), (D, B). Based on this ordered set, we can obtain two coefficients vectors c\, 
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Fig. 2 A dendritic structure 
used in the example calculation 
for the matrix method 




one for trips that begin at x and move towards B, denoted cj = (1 0000 0) r , and 
another for trips moving from x towards node A, denoted c\ = (0 1 0 0 0 0) r . Using 
the rules described in (16), we can construct the matrix Q for our dendritic structure 
as follows: 



Q = 



( o 10 
0 0 

0 0 

0 1 

0 0 

0 0 



_ 1 

3 

2 
3 

0 

2 
3 

V o 



3 

1 

3 

0 

7; 



Note that all rows and columns sum to 1 . Knowing this matrix Q and breaking the 
trips into four main classes, it is straightforward to find the complete Green's func- 



tion: 



G sg (x, y, t) = |^(e*~M) g G oo(A*(2n + x + y), t) 

+ J2(Q k ~ lc i) g G °°( Ax ( 2n +2-x+y),t) 
k=l 

^max 

+ J2(Q k ~ lc \)g> G °°( Ax ( 2n +2 + x-y),t) 
k=l 

^max 

+ J2(Q k ~ lc i) g >Goo(Ax(2n +4-x-y),t). 



(18) 



4 Convergence of Methods 

We first validate the computational implementations of our algorithms by construct- 
ing the Green's function Gij(x,y,t) on a small binary tree. Two profiles of the 
response function obtained by the length-priority, the Monte-Carlo and the matrix 
methods are shown in Fig. 3 and compared to a numerical simulation computed by 
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Time 

Fig. 3 The Green's function constructed by a number of methods. The voltage traces show the analyti- 
cal solutions Gij(x, y, t) for fixed x and y on a binary tree for the length-priority method (red circles), 
the Monte-Carlo method (blue diamonds) and the matrix method with & max = 100 (green squares) su- 
perimposed on NEURON'S numerical solution (black line). Parameter set A (Table 1) is used in these 
computations 

the software package NEURON [26]. These plots demonstrate an excellent agree- 
ment between a number of approaches for obtaining the Green's function and the 
numerical solution, with a slightly worse performance of the Monte-Carlo method 
for larger times. 

The computational convergence of any algorithm impacts both the accuracy of its 
results and the speed at which these results are obtained. While the series solution 
to the Green's function (5) is proven to converge for a sufficiently-high number of 
terms (see the Appendix), this assumes optimal ordering of terms in the solution. 
Because the coefficients A tr i p are impossible to compute until a trip is constructed, 
generating terms for the series solution in descending order of magnitude is inher- 
ently difficult. The four-classes, length-priority and matrix methods generate trips in 
order of increasing L tr i p , with the aim of ordering trips by their Goo(£trip, t) terms, 
which decreases monotonically in the length of the trip. The Monte-Carlo method 
uses a stochastic method to order trips by their probabilities, with more likely trips 
contributing more to (5). However, none of these approaches order the trips optimally, 
and hence their accuracy relies not on the theoretical convergence of the mathematical 
method, but on the computational convergence of the algorithm that implements it. 

To assess the computational convergence of the algorithms, Green's function so- 
lutions were constructed on a number of branching geometries shown in Fig. 4. In 
addition to a binary tree with the topology in Fig. 4A, the algorithms were also ap- 
plied to four real neuronal reconstructions obtained from the NeuroMorpho Database 
[27] in . swc format and shown in Figs. 4B-E. These files use a sequence of nodes 
with precise radii and three-dimensional locations to describe the location of the soma 
and the paths taken by the axon and each dendrite. A dendrite's path can be described 
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Fig. 4 Neuronal structures used in construction of the Green's function. A: a binary tree, B: a rabbit 
amacrine cell [28], C: a rat pyramidal cell [29], D: a rat Purkinje cell [30], and E: a blowfly tangential cell 
[31] 

using as many nodes as necessary to accurately reflect the spatial jitter and variation 
in radius of its path. Because radii are described at nodes, edges between two nodes of 
different radius taper. The sum-over-trips formalism requires constant diameter along 
edges, but allows discontinuous jumps in the diameters at nodes. Hence, edge diam- 
eter was defined as the average of the diameters of adjacent nodes. This allows full 
dendritic branches to be represented as a sequence of uniform cylinders of arbitrary 
length and with abrupt changes in diameters at nodes. 

We used the following normalised L 1 error as a measure of convergence: 

1 f 7 \ 

£ = tt \Gij(x,y,t)- V*(x,y,t) df, 
VN JO 

where T is the final simulation time, V*(x, y, t) is NEURON's numerical solution to 
very high accuracy and = fj V*(x,y,t)dt is the integral of the accurate NEU- 
RON solution. This convergence measure is therefore relative to the amplitude of the 
"real" solution, and thus errors e are comparable between different neuronal types. 

Figure 5 shows the convergence of the four-classes and of the length-priority meth- 
ods as a function of the number of trips on five geometries from Fig. 4. Three sets 
of parameters given in Table 1 were considered for the binary tree, and the relative 
errors s for each case are demonstrated in Figs. 5A-C. These plots illustrate fairly 
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Number of Trips Generated 



Fig. 5 Convergence of the four-classes and length-priority methods for a number of dendritic morpholo- 
gies. The relative error s of the approximation of Gij(x,y,t) is shown as a function of the number of 
trips in the sum-over- trips framework for injection at y and measurement at x on the dendritic trees in 
Fig. 4. Membrane parameters for real dendritic morphologies: C = 1 uFcm -2 , R = 3,000 £2 cm 2 and 
R a = 100 Qcm. Note that the four-classes method always begins with four trips, and each step in the 
algorithm adds a further trip of each class 

uniform convergence, in which both methods offer similar accuracies and rates of 
convergence. Of the two trees with biophysically realistic parameters, the binary tree 
with the longer branches (parameter set C) converges faster, as is expected on struc- 
tures with longer trips in each Class. This is reflected in Table 2, which shows the 
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Table 1 Parameter sets of the binary tree in Fig. 4 A 





Parameter set A 


Parameter set B 


Parameter set C 


Branch length L 


U.J 


50 jam 




100 um 




Branch diameter a 


0.05 


1 um 




1 um 




Diffusion coefficient D 


1 


2.5 x 10 4 


2 -1 
um ms 


2.5 x 10 4 


2 -1 

1 1 m ms 

JJ111 


Membrane time constant x 


1 


3.3 ms 




3.3 ms 




Membrane capacitance C 


1 


1 uFcm - 


2 


1 uFcm - 


2 


Table 2 Length-priority 
method on a binary tree: number 


Binary tree 


Relative error threshold e 




of trips required for a given 
accuracy 




0.1 


0.05 


0.01 


0.001 


Parameter set A 


3,240 


8,750 


1,820,000 


>5 x 10 7 




Parameter set B 


825 


2,600 


129,000 


>5 x 10 7 




Parameter set C 


22 


65 


815 


7,700 



number of trips required on the binary tree to remain under a given error threshold 
for different parameter sets. Moreover, a binary tree with non-dimensionalised pa- 
rameters (parameter set A) requires noticeably many more trips for desired accuracy 
in comparison to the same tree with biophysically realistic parameters. 

Figures 5D-G show s for the structures in Figs. 4B-E respectively. They demon- 
strate that convergence is non-trivial on complex branching structures. Figure 5D 
shows that the length-priority method makes consistently less error on the amacrine 
cell geometry, in contrast to the convergence of the Purkinje cell, shown in Fig. 5E, 
where the four-classes method generates less error for all numbers of trips. Both of 
these show strongly irregular convergence and high- amplitude oscillation in the er- 
rors 8 in the amacrine cell. For both methods, the Purkinje cell shows a plateau in 
error for Green's functions with few trips, indicating that either these trips are of 
small magnitude or that their voltage traces alternate between undershooting or over- 
shooting the correct solution between subsequent trips. This indicates that neither the 
length-priority or the four-classes methods are good heuristics for ordering terms in 
the Green's function. This is further hinted at by the oscillating property of the er- 
ror, which implies that there are regions where trips that increase the error are more 
frequent than trips that reduce it. 

The pyramidal cell's convergence shows very discontinuous behaviour (Fig. 5F), 
particularly in the length-priority method. The large jump in error when approxi- 
mately 350 trips are included in the Green's function was found to be caused by 
the first and shortest Class 2 trip included thus far, with all prior trips belonging to 
Class 1 . This behaviour is likely to arise if there exist very short branches along the 
shortest and most direct x -> y trip, and thus many Class 1 trips are generated first, 
being shorter than the first Class 2 trip. Whilst one of the motivating reasons for con- 
sidering a length-priority approach was to generate trips fully by length order, this 
heuristic makes no attempt to include the coefficient A tr ip in its ordering. This is an 
example of a pathologically large change in the coefficients value for a Class 2 trip 



4^ Springer 



Journal of Mathematical Neuroscience (2012) 2:11 



Page 19 of 28 



which contributes a very significant amount to the Green's function. The four-classes 
approach, which enforces generation of trips of all four classes at every added ex- 
cursion, does not show such a drastic drop in error. However, the error plot is still 
very discontinuous, and this may be a characteristic of situations as we have just de- 
scribed, where points x and y are placed on branches having a very different length 
to those on the most direct x —> y trip, or when these points are placed very close to 
a node. Whether injection and measurement points are located on branches that are 
significantly longer or shorter than those along the shortest x -> y trip, both the four- 
classes and the length-priority methods will generate trips in an "unnatural" order, 
subsampling the trips where current will spread the most, but oversampling in areas 
of the tree with very short branches. This pathological feature may not be inherently 
present in the real neuronal morphology, but may have been created during digital 
reconstruction from slice image data if, for example, a change of radius were found 
along the branch. Therefore, this pathology may not be representative of the neuronal 
geometry, but becomes a function of the reconstruction. 

The tangential cell's convergence, shown in Fig. 5G, shows almost identical errors 
for both the four-classes and the length-priority methods, indicating that trips are 
generated in a similar order regardless of method. Contrary to the example with the 
pyramidal cell, this behaviour is likely to occur when x and y are placed on branches 
that are significantly shorter than those that arise on the shortest x —> y path, such that 
the length-priority method returns trips of Class 1, 2, 3 and 4 in sequential order, as 
these increases in length are shorter than adding an excursion along the direct x —> y 
trip. 

Our results clearly indicate that the convergence of the realisation of the sum- 
over-trips framework by either the four-classes or the length-priority method strongly 
depends on a dendritic geometry. For real morphologies, the number of trips required 
quickly becomes very large to the point where guaranteeing convergence to within 
some small error threshold may become computationally expensive. 

The convergence of the Monte-Carlo method is shown in Fig. 6 for the binary tree 
in Fig. 4 A with parameter set A and for a larger binary tree of depth 16 with the same 
parameters. Boundary effects can be seen on the smaller tree, where the convergence 
rate is slightly faster than that of a typical Monte-Carlo integration observed here for a 
larger tree. It is worth noting that the x-axis on this plot shows the number of random 
walks generated; however, due to the number of sub-trips extracted from each random 
walk, the number of terms contributing to the Green's function can potentially be 
significantly different. The graphs show that the Monte-Carlo method is very slow to 
converge, although the method is much more predictable in its convergence despite 
the noise. As expected, therefore, the Monte-Carlo method generates trips that are 
more "naturally" ordered, and hence convergence is much more monotonic. Despite 
this improved ordering of terms in the series solution, the Monte-Carlo approach 
remains computationally intensive and very slow to converge with increased number 
of trips. 

Finally, the convergence of the matrix method on a binary tree is shown in Fig. 7. 
This algorithm converges extremely quickly to within very small error tolerances as 
a function of k max , the maximum number of edges covered by trips generated. The 
values of the product of A tr i p coefficients obtained by this method remain (9(1) for all 
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Fig. 6 Convergence of the Monte-Carlo method for a binary tree. The relative error s is shown as a 
function of the number of random walk realisations generated, k. A: error generated on the binary tree 
in Fig. 4A with parameter set A. The red line shows a fit for s ~ B: the convergence error on a 

binary tree of depth 16 (65,536 nodes). The red line demonstrates a fit for s ~ , the typical rate of 
convergence of a Monte-Carlo integration 




0 10 20 30 40 50 60 70 80 90 

/Cmax 



Fig. 7 Convergence of the matrix method. The relative error s is shown as a function of the maximal 
number of edges travelled in the trips, k max , for a binary tree in Fig. 4 A 

&max which agrees with the proof in [25]. Because the algorithm is based on simple 
matrix- vector multiplication, where the matrix is \£\ x \£\ in size, the computation 
of the Green's function for small trees such as the binary trees used here to within 
s = 10 -15 only takes a fraction of a second. On more complex trees, such as Purkinje 
cells, this becomes more expensive, although computing Gij(x, y, t) for the whole 
tree remains computationally preferable to the use of brute-force simulators. Using 
the reciprocal rule (9), this is possible in \6\ applications of the algorithm. This com- 
pares favourably with the length-priority and the four-classes methods, which require 
+ l)/2 applications of the algorithm, and with NEURON, which would re- 
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quire |£| 2 simulations, and this would only provide solutions for a single point y on 
each edge. In addition, the sparseness of the matrix Q means that coefficient calcu- 
lation up to & max only takes 0(\E\k max ) time, and so the method scales linearly with 
the number of branches on the tree. 

4.1 Structural-Electrotonic Properties 

The Green's function Gtj(x,y,t) constructed for a given dendritic geometry provides 
a measure of the transfer impedance between the input location y on branch j and 
the point of measurement x on branch i . Assessing whether the neuronal geometry 
significantly impacts this input-output relation is a step towards answering questions 
regarding the structure-function relationship behind the enormous natural variation 
in dendritic morphologies. Using the measures introduced by Zador et al. [32], the 
propagation delay and the log-attenuation, we analyse the transfer of the response 
signal in four reconstructed cells in Fig. 4. For a pair of points (x, y) along the tree, 
we reintroduce the propagation delay V xy as a measure of the impact of the tree's 
electrotonic structure on the timing of signals, defined as 

T^xy = tx ~ ty • 

Here, t x and t y are the centroids of the two corresponding voltage transients G x (t) = 
Gij(x, y, t) and G y (t) = Gjj(y, y, t) respectively: 

. _ f™tG x (t)dt - _ f™tG y (t)dt 

tx ~ f~G x (t)dt ^ h ~ f~G y (t)dt' 

This delay measure admits an additive property such that V xy = V xz + V zy for a 
point z between x and y. The log-attenuation of the response signal between a pair 
of points (x, y) is computed as C xy = \ogA xy > 0, where 

^ = { Q oo^T > 1- (19) 



f™Gy(t)dt 

f™G x (t)dt 



It acts as a measure of the amount a transient signal's amplitude diminishes as it 
travels between two points. C xy is also additive for a point z between x and y, that 

is? £xy = &xz ~l~ ^zy- 

Figure 8 demonstrates the propagation delay and the log- attenuation as a func- 
tion of distance y away from x for the reconstructed dendritic morphologies in 
Figs. 4B-E. The point x was placed near the soma as shown in Fig. 4 and the po- 
sition of y was moved away from x to the distal dendrites along a single path. As 
expected with an additive property, both the delay and the log-attenuation are lin- 
ear in the distance between y and x. The curves in Fig. 8 show a noisy linear trend, 
which could be smoothed to better demonstrate this linearity by sampling the data 
from multiple points located along different branches, but at the same fixed distance 
away from x, in the manner of an expanding sphere of radius y, with origin at x. To 
compare how the response signal is transferred in four neuronal types, we plot the rate 
of change of the delay and log- attenuation for individual cells in Fig. 9, computed by 
a linear regression and imposing that the line passes through the origin. It succinctly 
illustrates that the input signal will spread differently in these four cells, with the tan- 
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Distance (um) 



Fig. 8 Propagation delay and log-attenuation for reconstructed geometries 



Fig. 9 Transfer properties of 
the response signal for 
reconstructed geometries 
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gential cell having a similar rate of delay with the pyramidal cell and a similar rate 
of log-attenuation with the amacrine cell. The signal in the Purkinje cell is shown to 
be attenuated most, whereas the pyramidal cell transfers it very effectively. Similar 
conclusions, but about the propagation of the dendritic action potential, were made 
by Vetter et al. [30]. The activation of voltage-gated channels is expected to be less 
robust in the case of strong attenuation of the passive spread of voltage which might 
explain the results in [30] . 



5 Discussion 



In this paper, we introduced a number of efficient algorithms for the computational re- 
alisation of the sum-over-trips framework and assessed their convergence. We started 
with some modifications of the four-classes algorithm of Cao and Abbott [21] to 
avoid constructing duplicate trips. An unambiguous context-free grammar was de- 
rived, which is able to generate all trips uniquely and in monotonic order of length. 
We then developed the length-priority method, in which trips are constructed purely 
in length order rather than in classes. Both methods were found to demonstrate very 
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nonuniform convergence which was highly-dependent on the dendritic morphology, 
as well as the biophysical properties of the cell membrane. Oscillations of the con- 
vergence error make it difficult to predict the number of trips required for construct- 
ing the Green's function on a particular geometry. Dendritic structures with longer 
branches of uniform diameters will converge faster, that is, for a smaller number of 
trips in the series solution. Instead of sampling the trips in some well-defined order, 
we also derived a stochastic method of sampling the trips based on a Monte-Carlo ap- 
proach. Finally, we proposed an extremely efficient matrix method which computes 
the trip coefficients, Atrip, for trees where all branches are integer-multiples in length 
to some base length, Ax . 

Although we considered dendrites to be passive in this study, the proposed al- 
gorithms can be easily generalised to support quasi-active (resonant) dendrites with 
a calculation of the Green's function in the Laplace domain [24]. Moreover, it is 
straightforward to include an isopotential soma in the Laplace-domain series solu- 
tion. A soma can be considered as a special node with the factors pk(co) on the 

branches connected to the soma defined as pk(co) = r^ l J(co-\-z~ l )D^ 1 /(Cco + 



R -1 + J2 m r m 1 y + T ~ 1 )Dm )» where C and R are the capacitance and the resis- 
tance of the somatic membrane and co is the Laplace transform's frequency variable. 
Similar factors can also be found for when the leaky-end boundary condition, referred 
to as natural termination by Tuckwell [33], is imposed at the terminals. A knowledge 
of the Green's function for a given dendritic structure allows one to efficiently find 
the sub-threshold voltage response along the entire tree for any number of various 
inputs, either analytically or via a computation of the convolution integral. This obvi- 
ates the need for the brute-force numerical simulations of an underlining set of PDEs. 
Such simulations may be computationally expensive, particularly since they have to 
be re-initiated each time a new stimulus is introduced. In the case of supra-threshold 
inputs, which can activate voltage-gated channels known to be present in dendrites of 
many neurons, the Spike-Diffuse-Spike (SDS) type model [34, 35] can be utilised for 
analysing the propagations of dendritic action potentials. Although the voltage-gated 
channels in the SDS framework are modelled by piecewise linear instead of nonlin- 
ear dynamics, it has been shown that the speed of a wave propagation in the SDS 
model is in excellent agreement with a more biophysically realistic nonlinear model 
[36]. However, an analytically tractable SDS model combined with a fast algorithm 
for constructing the Green's function on real geometries provides a computationally 
efficient framework for studying wave scattering in dendrites. 

Although networks of spatially-extended neural cells can be numerically simu- 
lated, there are currently few mathematical studies of such networks. A natural exten- 
sion might be to consider a network of branched neurons coupled by gap-junctions. 
The sum-over-trips formalism can then be generalised to support a presence of new 
boundary conditions. Recent results of Harris and Timofeeva [37] can be applied to 
the case of tip-to-tip coupling of the dendritic branches. The proposed algorithms 
can then be modified by including additional sum-over- trips rules. It is worth men- 
tioning that the computational schemes presented in this paper are able to handle 
cyclic graphs, which may form as a result of gap-junction coupling across several 
neurons. While the matrix method is expected to be the most efficient for a network 
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of symmetric or regular structures, realistic reconstructions of discretised trees re- 
main within computational reach. For example, given a sparse random matrix, of 
correct density and of size 20,000, equivalent to a dendritic tree with 10,000 edges, 
the calculation of Gij(x,y,t) for all j up to & max = 10 5 only takes fifteen seconds 
on a desktop computer. For comparison, the Purkinje cell reconstruction in Fig. 4D 
has just under 5,000 branches. For the case of very large, complex irregular struc- 
tures it might be possible to employ a recently developed technique of reducing the 
complexity of large dendrites [38] before applying the sum-over-trips methodology. 
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Appendix: Mathematical Convergence of the Sum-over- trips Series Solution 

Here, we consider an identical diffusion coefficient D for all branches, although it is 
possible to generalise this proof to support different diffusion coefficients. Fixing t 
throughout, we let 

oo 

Gij (X, y) = ^tripGoo(Xtrip) = ^ ^ ^ trip ^oo(^ trip), 

trips k=0 paths with 

k nodes 

where 

GooCLMp) = -^e-^e^. (20) 

A tr ip is a product of k factors 2p e (0, 2) and 2p — 1 E (- 1,1), where k is the number 
of nodes visited by the trip, for any branch diameters. Then, for a trip touching k 
nodes, 

|A t rip|<2*. (21) 
There exists a constant B > 0 such that every trip touching k nodes satisfies 

L trip > Bk. (22) 
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This makes B the coefficient of the lower bound on trip length in terms of the number 
of nodes in a trip. Intuitively, Bk equals the minimum distance between any two 
nodes where, for this purpose, we count x and y as nodes. 

Let Fk be the number of trips with k nodes. Since each node has degree d < 3, 
then 



We introduce 



Then 



F k <3 k . 



He — ^trip G(x> (X trip)- 



trips with 
k nodes 



\r k \ 



^ * ^ trip Goo (X trip) 



trips with 
k nodes 



< 



^ ^ |^tripGoo(Lt r ip) | 



trips with 
k nodes 



^1 l^tripl |Goo(L tr ip)|. 



trips with 
k nodes 



(23) 



(24) 



—EL 

For simplicity, we rewrite the Green's function (20) as (L trip ) = Ce tri P, where 

1 



c = 



Then using (21)-(24) we get 



; and E = 



\r k \< 2 * Ce 



-ELt 



trips with 
k nodes 

< £ 2*Ce" 

trips with 
k nodes 

< F k 2 k Ce~ EB2k2 

< 3 k 2 k Ce~ EBlk2 
= 6 k Ce- EB2k2 
= Ce- 



il B-k 



trip 



2 1-2 



(25) 
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We define N = |_ln(6) / (E B 2 ) J such that EB 2 k - ln(6) > 0 for Wk > N. Then 



£=0 



(X) 

£=0 

oo 

< \^ Ce -/:(£5 2 /:-ln(6)) 
iV 



£=0 £=N+1 



(26) 



The first sum in (26) is a finite sum of finite terms, and is hence finite. We will 
now show that the second sum is also finite using d'Alembert's ratio criterion for 
convergent series. The ratio pk of the consecutive terms in the series, k and k + 1, is 



Pk = 



Ce -(£+l)(£5 2 (£+l)-ln(6)) 



C Q -k{EB 2 k-\n{6)) 



= e -(£5 2 (2^+l)-ln(6)) > 



Letting k -> oo, we obtain 



Poo = lim p k 

k^oo 



. \{ m e -(£5 2 (2/:+l)-ln(6)) 
k^oo 



= 0. 



With 

Poo < 1? the second sum in (26) converges absolutely for all constants 
B, C, E > 0. Therefore, the series in (26) is absolutely convergent for sufficiently- 
high k. 

If we define 

M 

GfJ(x,y,t) = J2 A tripGoo(Ltrip,0, 
k=0 trips with 
k nodes 

then 

oo 

\Gij-Gff\< Ce-W* 2 *- 1 *™ 

k=M+\ 

and the path integral converges faster than e - ^ in the worst case, with the number of 
nodes k visited by the trips. 
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